Mars Declination Periodogram

This notebook plots the declination of mars, and transforms it into the frequency domain.

We use the PyEphem module to calculate mars positions.

Animation depends on the external utility ImageMagick



In [1]:

    
import ephem
import pandas as pd

Build up a dataset of mars posiitons



In [3]:

    
mars = ephem.Mars()

#these are kluge functions to get the data in the format we want. Not efficient!
def getdec(date):
    mars.compute(date)
    return mars.dec

def gethlong(date):
    mars.compute(date)
    return mars.hlong

def getsundistance(date):
    mars.compute(date)
    return mars.sun_distance

mars_data = pd.DataFrame(index=pd.DatetimeIndex(start='Jan 1, 1700', end='Dec 31, 2260', freq='M'))
mars_data['date'] = mars_data.index
mars_data['decl'] = 180./pi*mars_data.date.apply(getdec)
mars_data['hlong'] = mars_data.date.apply(gethlong)
mars_data['sun distance'] = mars_data.date.apply(getsundistance)
mars_data['year_fraction'] = mars_data.index.year + 1.*mars_data.index.dayofyear/365.
mars_data.head(1)









    Out[3]:






  
    
      
      date
      decl
      hlong
      sun distance
      year_fraction
    
  
  
    
      1700-01-31
      1700-01-31
      -12.707192
       3.168863
       1.641473
       1700.084932

Plot the declination

We'll also draw (by eye) a curve that follows the location of the declination event in a 15 year period.



In [4]:

    
#eyeball fit these parameters
a = 15.75
b = 10.4
mars_data['retroprog'] = mars_data['year_fraction'].apply(lambda x: .4*180./pi*sin(2.*pi*(x+b)/a))

plt.figure(figsize=(10,3))
mars_data['decl'].loc['1964':'2014'].plot(label='Mars Declination')
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.xlabel('Year', fontsize=18)
plt.ylabel('Declination [Deg]', fontsize=18)
plt.title('Declination of Mars', fontsize=18)
mars_data['retroprog'].loc['1964':'2014'].plot(color='c', label='Retrograde Progression', linewidth=3, alpha=.5)
plt.legend()









    Out[4]:





<matplotlib.legend.Legend at 0x10b9bca50>

Periodogram

Now lets calculate the power spectral density, and draw a periodogram.

We transform the data from units of frequency to units of period.



In [5]:

    
import scipy.signal

freq, power = scipy.signal.periodogram(mars_data['decl'], fs=12, scaling='spectrum')

plt.figure(figsize=(10,3))
plt.plot(1/freq, power, '-', linewidth=5, alpha=.5)

plt.xlabel('Period [Years]', fontsize=18)
plt.ylabel('Power Spectrum Density', fontsize=18)
plt.ylim(0, 200)
_, _, ymin, ymax = plt.axis()
plt.vlines(x=1, ymin=ymin, ymax=ymax, colors='g', label='Earth Orbital Period')
plt.vlines(x=1.8808, ymin=ymin, ymax=ymax, colors='r', label='Mars Orbital Period')
plt.vlines(x=15.75, ymin=ymin, ymax=ymax, colors='c', linewidth=3, 
           alpha=.5, label='Alignment of Retrograde with Mars Perigee')
plt.xticks([1,1.8808, 5, 10, 15.75, 20],[1,1.88, 5, 10, 15.75, 20], fontsize=14)
plt.yticks(fontsize=14)
plt.title('Power Spectrum Density of Mars Declination Measurement', fontsize=18)
plt.legend(fontsize=14);
plt.ylim(ymin, ymax)
plt.xlim(0,18)









    Out[5]:





(0, 18)

Mars phase plot

Now we'll do a coordinate transform to isoate the x and y axes



In [7]:

    
mars_data['x'] = mars_data.apply(lambda row: row['sun distance']*cos(row['hlong']), axis=1)
mars_data['y'] = mars_data.apply(lambda row: row['sun distance']*sin(row['hlong']), axis=1)
mars_data['dx'] = np.append(np.diff(mars_data['x']), 0)
mars_data['dy'] = np.append(np.diff(mars_data['y']), 0)
mars_data.head(1)









    Out[7]:






  
    
      
      date
      decl
      hlong
      sun distance
      year_fraction
      retroprog
      x
      y
      dx
      dy
    
  
  
    
      1700-01-31
      1700-01-31
      -12.707192
       3.168863
       1.641473
       1700.084932
      -13.728108
      -1.640862
      -0.044758
       0.071736
      -0.355585

Generate the phase plot data



In [8]:

    
spacing = .2 #how far apart should the quivers be.
x, y, dx, dy = np.meshgrid(np.arange(-2.25, 1.75,1.25*spacing), np.arange(-1.75,2.25,1.25*spacing), 
                           np.arange(-.75, .75, spacing), np.arange(-.75, .75, spacing))

a = -1.0
r2 = x**2 + y**2
theta = np.arctan2(y,x)
ddx = a/r2*np.cos(theta)
ddy = a/r2*np.sin(theta)

Draw a non-animated plot.

This actually sets the x and y lims for the animated plot, becuase those are harder to calculate when you are using a subset of the data. Its a bit of a sloppy coding idea, but hey...



In [10]:

    
fig = plt.figure(figsize=(12,12))
#rc('text', usetex=True)

ax1 = plt.subplot2grid(shape=(3,3), loc=[1,0], rowspan=2, colspan=2)
plt.plot(mars_data['x'], mars_data['y'], alpha=1)
plt.plot(0,0, 'yo', markersize=20)
plt.axis('equal')
plt.xlabel('X Position  $\propto$ X Potential Energy', fontsize=16)
plt.ylabel('Y Position  $\propto$ Y Potential Energy', fontsize=16)
xmin, xmax, ymin, ymax = plt.axis()
plt.box('off')
#plt.title('Physical Position', fontsize=18)


####### x vs x velocity
plt.subplot2grid(shape=(3,3), loc=[0,0], colspan=2, rowspan=1)

horiz_pos = x[0,:,:,0]
vert_pos = dx[0,:,:,0]
horiz_change = ddx[0,:,:,0]
vert_change = dx[0,:,:,0]

plt.quiver(horiz_pos, vert_pos, vert_change, horiz_change, scale=6, headlength=10, headwidth=6)

plt.plot(mars_data['x'][:-1], mars_data['dx'][:-1], alpha=1)
plt.ylabel('X Velocity  $\propto \sqrt{X\,Kinetic\,Energy}$', fontsize=16)
#plt.axis('equal')
plt.box('off')
plt.xticks([])
plt.title('X position vs X velocity', fontsize=18)
plt.xlim(xmin, xmax)
plt.ylim(-.5, .6)

####### y vs y velocity
plt.subplot2grid(shape=(3,3), loc=[1,2], rowspan=2, colspan=1)

vert_pos = y[:,0,0,:]
horiz_pos = dy[:,0,0,:]
vert_change = ddy[:,0,0,:]
horiz_change = dy[:,0,0,:]

plt.quiver(horiz_pos, vert_pos, vert_change, horiz_change, scale=2, headlength=20, headwidth=12)

plt.plot(mars_data['dy'][:-1], mars_data['y'][:-1], alpha=1)
plt.xlabel('Y Velocity  $\propto \sqrt{Y\,Kinetic\,Energy}$', fontsize=16)
#plt.axis('equal')
plt.yticks([])
plt.box('off')
plt.title('Y position vs Y velocity', fontsize=18)
plt.ylim(ymin,ymax)
plt.xlim(-.5, .5)


fig.subplots_adjust(wspace=0, hspace=0)

Now lets animate it



In [12]:

    
from matplotlib import animation

def animate(nframe):
    nframe *= 1
    fig = plt.gcf()
    plt.cla()

    ax1 = plt.subplot2grid(shape=(3,3), loc=[1,0], rowspan=2, colspan=2)
    plt.plot(mars_data['x'][:nframe], mars_data['y'][:nframe], alpha=1)
    plt.plot(mars_data['x'][nframe], mars_data['y'][nframe], 'ro', markersize=10)
    plt.plot(0,0, 'yo', markersize=20)
    plt.axis('equal')
    plt.xlabel('X Position', fontsize=16)
    plt.ylabel('Y Position', fontsize=16)
    #xmin, xmax, ymin, ymax = plt.axis() #set these in the one-off run
    plt.xlim(xmin, xmax)
    plt.ylim(ymin, ymax)
    plt.box('off')
    #plt.title('Physical Position', fontsize=18)


    ####### x vs x velocity
    plt.subplot2grid(shape=(3,3), loc=[0,0], colspan=2, rowspan=1)

    horiz_pos = x[0,:,:,0]
    vert_pos = dx[0,:,:,0]
    horiz_change = ddx[0,:,:,0]
    vert_change = dx[0,:,:,0]

    plt.quiver(horiz_pos, vert_pos, vert_change, horiz_change, scale=5, headlength=10, headwidth=6)

    plt.plot(mars_data['x'][:nframe], mars_data['dx'][:nframe], alpha=1)
    plt.plot(mars_data['x'][nframe], mars_data['dx'][nframe], 'ro', markersize=10)
    plt.ylabel('X Velocity', fontsize=16)
    #plt.axis('equal')
    plt.box('off')
    plt.xticks([])
    plt.title('X position vs X velocity', fontsize=18)
    plt.xlim(xmin, xmax)
    plt.ylim(-.5, .5)

    ####### y vs y velocity
    plt.subplot2grid(shape=(3,3), loc=[1,2], rowspan=2, colspan=1)

    vert_pos = y[:,0,0,:]
    horiz_pos = dy[:,0,0,:]
    vert_change = ddy[:,0,0,:]
    horiz_change = dy[:,0,0,:]

    plt.quiver(horiz_pos, vert_pos, vert_change, horiz_change, scale=1.5, headlength=20, headwidth=12)

    plt.plot(mars_data['dy'][:nframe], mars_data['y'][:nframe], alpha=1)
    plt.plot(mars_data['dy'][nframe], mars_data['y'][nframe], 'ro', markersize=10)
    plt.xlabel('Y Velocity', fontsize=16)
    #plt.axis('equal')
    plt.yticks([])
    plt.box('off')
    plt.title('Y position vs Y velocity', fontsize=18)
    plt.ylim(ymin,ymax)
    plt.xlim(-.5, .5)

    fig.subplots_adjust(wspace=0, hspace=0) 

    
fig = plt.figure(figsize=(12,12))

anim = animation.FuncAnimation(fig, animate, frames=100)
anim.save('demoanimation.gif', writer='imagemagick', fps=5);